function x_new = rk4_resistive(x0, dt, b, r, c, n, lambda_R, I_Chaos, VPC)
    % rk4_resistive - 使用四阶 Runge-Kutta 方法对阻抗耦合模型进行单步更新
    %
    % 输入:
    %   y0       - 当前状态向量 [V1; iL1; V2; iL2] (应为列向量)
    %   dt       - 单步步长
    %   b        - 无量纲参数，通常为 E_j/V0
    %   r        - 无量纲参数，通常为 R_j/a
    %   c        - 无量纲参数，通常为 a^2 * C_j/L_j
    %   n        - 与压电传感器相关的参数
    %   lambda_R - 阻抗耦合强度（lambda_R），如论文中 lambda_R=10
    %   I_Chaos  - 混沌电流（无量纲）
    %   VPC      - 外部声信号（压电电压）
    %
    % 输出:
    %   y_new    - 单步更新后的状态向量 [V1; iL1; V2; iL2]
    
     x0 = x0(:); % 强制转为列向量

    k1 = resistive(x0, b, r, c, n, lambda_R, I_Chaos, VPC);                         % 第一阶段（k1）
    k2 = resistive(x0 + 0.5*k1, b, r, c, n, lambda_R, I_Chaos, VPC);                % 第二阶段（k2）
    k3 = resistive(x0 + 0.5*k2, b, r, c, n, lambda_R, I_Chaos, VPC);                % 第三阶段（k3）
    k4 = resistive(x0 + k3, b, r, c, n, lambda_R, I_Chaos, VPC);                    % 第四阶段（k4）

    % 更新状态
    x_new = x0 + dt/6*(k1 + 2*k2 + 2*k3 + k4);       % 使用RK4公式更新状态
end


function dXdt = resistive(X, b, r, c, n, lambda_R, I_Chaos, u_PC)
    % resistive - 计算阻抗耦合模型的状态导数
    % 输入:
    %   X      - 当前状态向量 [V1; iL1; V2; iL2]
    %   b, r, c - 无量纲参数，分别对应 E_j/V0, R_j/a, a² C_j/L_j
    %   n      - 与压电传感器相关的参数
    %   lambda_R - 阻抗耦合强度
    %   I_Chaos - 混沌电流（无量纲）
    %   u_PC   - 外部声信号（压电电压）
    %
    % 输出:
    %   dXdt   - 状态导数 [dx1; dy1; dx2; dy2]

    x1 = X(1);
    y1 = X(2);
    x2 = X(3);
    y2 = X(4);

    % 根据论文方程（7）计算各变量的导数：
    dx1 = x1 - (1/3) * x1^3 - y1 + I_Chaos + lambda_R * (x2 - x1);
    dy1 = c * (x1 + b - r * y1);
    dx2 = x2 * (1 - n) - (1/3) * x2^3 - y2 + n * u_PC - lambda_R * (x2 - x1);
    dy2 = c * (x2 + b - r * y2);

    dXdt = [dx1; dy1; dx2; dy2];
end

